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I. INTRODUCTION 


Two important parameters of interest in the field of aero- 
dynamics of airfoils are lift and drag. These can be evalu- 
ated either experimentally or theoretically. The desire for 
computational methods to aid the design process is promoted 
by the need to reduce the number and cost of wind tunnel tests. 

This thesis treats the problem of incompressible, two- 
dimensional steady flow about airfoils or airfoil combinations 
at large angles of attack. Such flows exhibit strong viscous 
flow effects including regions of flow separation. Therefore 
methods are required which can account for these effects. 

Currently there exist two main methods, namely the direct 
computation of viscous flows by means of the Navier-Stokes 
equations or the so-called viscous/inviscid interaction method. 
The former approach is more straightforward but also much more 
expensive and time-consuming. Therefore, the latter approach 
is to be preferred if it can be shown that it produces good 
agreement with the available experimental results. 

It is the purpose of this thesis to contribute to the 
evaluation of the viscous/inviscid interaction method. To 
this end, the viscous/inviscid computer codes developed by 
Cebeci and collaboratorsat the Douglas Aircraft Company were 
obtained and applied to several airfoils. 

In addition, a separate panel method was formulated and pro- 
grammed in order to obtain the inviscid flow over airfoil 
combinations. 
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The basic equations are formulated in Chapter II. Chapter 
III addresses the problem of inviscid flow calculations using 
the panel method. In Chapter IV the solution of the boundary 
layer equations by means of the Cebeci-Keller box method is 
explained. Finally, Chapter V describes the viscous/inviscid 


interaction problem and presents results of computations. 
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ЇЇ. BASIC EGUATIONS 


A. INTRODUCTION 


In this chapter, the basic equations of fluid flow are 
derived. We find that the resulting equations are PDE's whose 
exact solutions exist only in very few cases. The PDE's are 
classified, “parabolic, miny perDo e Uc E ела 
on the geometry of their domains of dependence and regions of 
influence, and the solution procedures are different according 


to the type of equation. Table 2-1 gives a brief classification 


of these equations. 


TABLE 2-1 


GLASSIFPICATION OF РОБ 


| EG SEM Parabolic | Hyperbolic 
| | 
J 


| 





Physical Upstream | No Upstream No Upstream 

Meaning Influence | Influence Influence 

гхапрїе -Laplace e Thin Shear . Supersonic 
Equations Layer Flow 


| 
- Steady Navier- 
Stokes 


Бие perturDation poiat 





Hx is domain on which solution 


at P depends 


b region of influence of 
а гепапкпасвстот ае в 


For example, the Laplace equation is elliptic. Its solution 
would have to be repeated for many iterations so that the up- 
stream influence can be gradually propagated (panel method in 
Chapter III), but the thin shear layer equations are parabolic. 
THicimenumer teal solution can be obtained by marching in the 


downstream direction only (Box method in Chapter IV). 


В. DERIVATION OF GENERAL EQUATIONS 

The continuity equation and the Navier-Stokes equations 
are basic for an aerodynamic analysis. We start with the basic 
physical concepts and derive the general equations for 2-D, 
unsteady, compressible, viscous flow. 

ime Continuity Eguation 

One of the basic laws of "Newtonian mechanics" states 
that mass can neither be created nor destroyed. Therefore, 
for a fixed control volume (see Figure 2.1), the principle of 
mass conservation can be stated that the net mass flow rate 
into and out of the control volume equals the time rate of 
change of mass within the control volume. 

If the central point 'P' has representative fluid proper- 
ties (velocity, density, pressure, etc.), then properties at 
other locations can be obtained by Taylor series expansions. 
Therefore the x-component of the velocity at the center of the 
positive x-face (right-hand face) is 


дїї Cl 


Du Ges 2 
Du 3x 2? "ied (2.1) 
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r5 





и + 


«e [o 
хіс 
ВЭ 

>< 


о г р u 


; | 


ал е 


Figures2.1.°8ConGrol Volunerrtons2—» 


As dx goes to zero, all higher order terms will be dropped, 
so that only the first two terms will be considered. Similarly; 


the density at the center of the positive x-face is 


Q2 


р dx, . (2.2) 


НЫ 22 


Then the mass flow rate out of the positive x-face is 


Mass Flow Rate Out (Density) (Velocity) (Area) 


др dx ах 


= (р 2 357: 2 ыа! ЫШ ЕЕ ДЕ mz а s, a cuin 
д ах 
= [pu + 7x (Pu) =I dy (2.3) 


By the same procedure the mass flow rate into the control volume 


through the negative x-face (left-hand side) is 
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Mass Flow Rate In = [pu - а. (ри) 3X] dy АА, 
From Eqs. (2.3) and (2.4), we get the net mass flow rate 


through the control volume in the x-direction. 
а к ах б Dr dx 
Net Mass Flow Rate = [pu yx (PU) I] dy [pu + yx (PU) =] dy 
== а 82:25) 
5 Х Y : ° 


In a Similar manner, the net mass flow rate in the y-direction 


is 
- 2 (ov) ax ду DN 
ду 


The total mass flow rate through the control volume is 


obtained by summing Eqs. (2.5) and (2.6). 
Total Net Mass Flow Rate = -[S (pu) ш (ov) 1x dy (2.7) 


Next, we consider the time rate of change of mass within 


the control volume. 


Time Rate of Change _ од 
of Mass c и 
- 92 ахау (2.8) 
dt 


У 


Now we combine Еав. (2/7) andi (2.3) using Ере сошес оно 


conservation of mass. Then 


үен Om _ Эр 
2000 cr уу (У) 14х dy = =) dxdy (2.9) 


Therefore we obtain the general form of the continuity 


сача топ for two-dimensional flow as 





д (рч) д (ру) Әр _ 
а + то + 5 = 0 (2.1503) 
For steady or unsteady incompressible flow, Eq. (2.10) reduces to 
Qu OV с 
D бу - 0. (2.10239) 
2. Navier-Stokes Equations 


Newton's second law of motion, when mass is conserved, 
equivalently states that the rate of change of the momentum of 


a body equals the sum of the forces applied to that body, or 


О, 


|] Р = 500) (2.11) 
In considering a small volume element of fluid, there are two 
types of forces to be considered, namely surface forces which 
are acting on the surface and body forces which are acting on 
the fluid inside the elemental volume, such as gravity (see 


PIgure 2.290 
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Figure 2:2. 
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Forces Acting on the Fluid 


Assuming that the stresses are known at point 'P', we get 


expressions for the stresses on the fluid element surfaces 


by Taylor series expansion. 


Net force in 
=d rection due 


to Normal Stresses 


Net force in 
x-direction due 
to Shear Stresses 


Therefore the total 


formed by summation 


И 
а 





|| 
а 





surface 


Of EGS. 














90 9 
Щас XK dx 
и ау 
ах ау (22:12) 
90 X G 9g ха 
ще к, к ах - [0 = аў СУ ах 
ду 2 ух oy 2 
dxidy 2213) 


forces in the x-direction are 


(2 ГО) М апс (213). 


С 90 
хх | ух 
OX ду 





) f, (surface) - | ] ахау (2.14) 


Let f (body) be defined as the body force per unit mass with 


the following components: 
= = : 
Г (body) =I EE 
Thus, the x-component of the body force is 


= (body) 


il 
3 
> 


о дхау.1]1.х (2. Е 


Adding Eqs. (2.14) апа (2.15) provides” the total force ined 


x-direction. 





) B = те (surface) and £ (body) 
9g 90 х 
= [ox + Era T ду 1 ахау (2:1 0) 


Now we consider the rate of change of the momentum of the fluid. 


Let us take the x-component only. Then, since the mass is constant 


d > _ йи 
ae TONE Ма 
_ du ди ди 
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ВЕНЕ, uU = (Е). УЕ), Е], апа Бу сһҺа1п-га1е 


du _ ды ах | да dy | 9u 9t 
dit с караны дуза WEE ot 
Е дъ дъ дъ 
B" 3x2 V ov "3t 


Substitution of Eqs. (2.16) and (2.17) into the x-component 


of Eq. (2.11) gives the final equation of motion. 


да 90 
XX | ух 
9х ду 





| y Een cou 


Don JX ay 7 9% 
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Now we want to show the stress in terms of the velocity com- 
ponents. In this thesis we will consider only simple "Newtonian" 
fluids obeying Stokes' law. This means that the 'extra' stress 
(above hydrostatic pressure) is proportional to the rate of 


strain. With the definition, 


Extra Stress = constant x (rate of strain) 
and introducing y = coefficient of viscosity, 
B Qu 
Охх 5 S 29 im 
ди QV 
= = == ---- 
oxy и Зу OX 


where the pressure in an incompressible fluid is seen to be 


equal to minus one-third the sum of the three normal-stress 
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components in view of Eq. (2.10a). In two dimensions then, 


О + с = -2Р. Eq. (2.18) then becomes, if the body force 
хх уу 


is neglected 





9ч ,449u ,4, 98e — iP AT xx , 1 ху (2.19) 
9 дх оос Жр р 


Substitution then produces, for incompressible flow, 


2 2 
du du du 1 ӘР ð u g u 
dt 9х ду о ðX T ase 
where v = п/р = kinematic viscosity, and similarly for the 
y Component 
д OV ду oP Be За 
низ чута = - 2 + У + ol (2.200) 
x y ж эх ду 


These are the well-known Navier-Stokes equations for two- 


dimensional incompressible viscous flow. 


C. INVISCID FLOW EQUATION 

All real fluid flows are viscous, but inviscid flow can 
be assumed outside of a thin boundary layer and a narrow wake 
behind the body for large Reynolds numbers. This is the reason 
why the inviscid flow equations are important even though they 
represent an ideal case. If the flow far upstream is uniform 
then it is also irrotational. This allows the introduction of 


the velocity potential. 
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1. Velocity Potential 


The velocity potential is a function whose gradient 


Peeesents the fluid velocity. Thus 


Qo 
E 
| 
о? 
сэг 


| 
г 
| 
| 


= у (22792727 


Qo 
= 


where 


Qu zu px E 


Piemerore, the Significance of the velocity potential lies 
in the fact that one equation for $ can be used rather than 
three equations for the velocity components. 
2. Laplace Equation 
For steady, incompressible flow, the continuity equa- 


tion (2.10) becomes 


DIES ОМ ы у з. (222: 


This equation can be written in terms of velocity 


КОЕ Па! > by substituting Eq. (2.22). Thus 
2 2 
Әна. оар. б (2.24) 
2 2 
х ду 
This is the well-known Laplace equation (vector form is un = i.) 


It is a linear equation which makes it possible to apply the 


principle of linear superposition. For instance, 
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If Фү, Pos axon on satisfy 18 = 0, then also 


ф 91 t Ф, то. $5 satisfies Уф = 0. 


Thus the flow resulting from the superposition of incompressi- 
ble, irrotational flows is also an incompressible and irro- 
tational flow. This superposition principle makes it possible 
to build up quite complicated flow from a few simple solutions 
Of Laplace's equation. The singularity (or panel) methods 


presented in the next chapter are based on this idea. 


D. THIN SHEAR LAYER EQUATIONS 

High Reynolds number flows over airfoils (and other con- 
figurations) generate a very thin shear layer (boundary layer) 
close to the airfoil surface. If ô denotes the boundary layer 
thickness and x the distance from the leading edge of a flat 


plate, then it is well-known that 





S (x) ~ Vvx/U ог 2-1 све 
x X 
where v is the kinematic viscosity and Re, = Ux/v. This formal 
shows that 
шонг NET if не аа. 
х x 


Hence the flow outside of the boundary layer can be considered 
to be inviscid, but the effect of viscosity cannot be neglected 


within this layer. Nevertheless, the Navier-Stokes equations 


24 


for steady incompressible flow can be simplified because of 
the fact that 6 is much smaller than the representative length 
шеЩешеталт вой (the chord length). This canbe sdeducedmfrom 


the Navier-Stokes equations as follows: 


2 2 
1 әр ð u ð u Qu ди 
- = — + у(-5- + 6) с а чу — (2725) 
р дх 3x^ o 9х ду 
2 2 
EET OP ду 9 Уу ” OV OV 
У 9х ду У 


u is now replaced by a typical value, say О 
X is now replaced by a typical value, say 2; 


y is now replaced by a typical value, say 6. 


Then 5з can be expressed by О /9; 
9 


u 
ax сап be expressed by Ue t: 


ОР 
эх 





can be expressed by pus /2 
e 


(because P and Ue are related by the Bernoulli equations). 
Therefore the x-component terms of the Navier-Stokes equation 


can be estimated to have the following magnitudes: 


2 2 
1 ӘР ð u ð u дъ Qu 
--д- + v (— + —) = u = t v >= 
р 9х эх? 72 эх ду 
2 2 2 
"e е е 
ди Г И 


25 


du 


The magnitude of the term v 5v follows from the application 
U 
of the continuity equation ou + 27 Опет _— ae Sor or 
8U Ox дү 1 9 
маш -- and therefore 
UJ i 
у 9% ЕЕ с 
ду Ll ô & 


The two viscous terms are of vastly different magnitude 


because 21114 


ч U /1* and 3 ^u/3y* Е u/s * hence 

2 2 2 2 2 2 

д ц/ду >> д ц/дх апа д u/dx can therefore be neglected 
compared to 3^u/3y^. Finally, the term v 3^u/3y^ must be of 
the same magnitude as the other terms if the influence of 
viscosity is to be retained. The v-compenentotermsmer te 


Navier-Stokes equations are easily estimated to be smaller 


than the x-component terms because 


92 
9х е 12 ду 02 
and hence are smaller by a factor 5. Therefore the two 
equations reduce to 
gu дц 1 dP < 
ал КУ; = - Se tv as (2.2% 
Y р y 
ОР _ 
ве = 0 (2.26 


By adding the continuity equation (Ес 2222 ксі гресе те 
tions, we get the basic equations to describe laminar flow 


thin shear layers. 
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E. TURBULENT FLOW 


We must deal with the instantaneous properties in the 


turbulent flow. Thus the time-mean value concept is applied: 


where u is the time-mean value, and u' is the fluctuating 


КОР ТОП. Similarly, 


Introducing these relations into Eq. (2.20) 
— - 2 2 
= ðu - (12 l OP д д 
Ше мая = - + (+ фа 
9х ду р 9х T 12 
ou ui да "у" 
ат (222229) 


We can see that pu'u' and pu'v' correspond to a normal stress 


and a shear stress. We call these terms turbulent stresses 


or Reynolds stresses. 
Similar analyses can be done for the y-component equations 
and z-component equations in the three-dimensional case. The 


extra Reynolds stresses can be summarized by the following 


array, 


21 




















О udi ub 
хх oxy O xz = 
e 2 

О @ g ccc dg. ML u'w' 

ух УУ yz TUNE 
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(Т С ! 1 uus um 

2Х ZY 922 mw 


Much of the effort in turbulent flow studies centers on the 


proper modelling of these turbulent terms. 
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ТТТ. РАПЦЕН МЕШЕСЕ 


А. INTRODUCTION 

The panel method was developed in the 1960's at McDonnell 
Douglas by Smith and Hess as a numerical approach for 2-D 
and 3-D potential flow problems. This chapter presents the 
application of the panel method to 2-D problems about one or 
two bodies. The basic idea consists in representing the flow 
past a body by a distribution of singularities (sources, 
sinks, vortices) on the surface such that the body surface 
becomes a streamline. 

The numerical approach requires some approximations (the 
assumption given in parentheses refers to our approach). 


A. The surface of the body is replaced by a finite 
number of elements (straight-line-elements). 


Е? tne condition of tangential flow is satisfied at a 
Ши эг number of points, the so-called control- 
points (midpoints of elements). 


С. The singularity distribution of each element is 
approximated by some kind of analytical functions 
(singularity strengths are assumed to be constant 
along any one element, but vary from element to 
element). 


The advantages of the panel method in comparison with other 
pasccagmres are: 


A. The panel method does not include an approximation in 
све physics--thin airfoil theory does. 


B. The panel method can be easily applied to both 2-D 
and 3-D problems--a virtually unsolvable task for 
conformal mapping procedures, which are confined 
to 2-0 configurations only. 
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C. The panel technique can be readily extended to flow 

fields past several bodies--a task which causes at 

least some troubles in "classical" mapping techniques. 
The method's versatility has been proven in various extensions, 
e.g., hydrofolls, cascades, nozzles, and even complete air- 
craft. Since its origin the method has been improved by 
features like higher order approximations to both body surface 
and singularity distribution, taking account for compressi- 
bility effects, and inclusion of wake models. Today the panel 


method is probably our most powerful tool in analyzing poten- 


tial flows. 


В. NONLIFTING FLOW PAST A BODY 

The effects of lift (respectively, camber and angle of 
attack) and displacement (resp. thickness) can be studied 
separately because of the linearity of Laplace's equation. 
This section is concerned about displacement flows due to the 
thickness of bodies, a flow which is usually represented by 
sources and sinks. 

We will first draw the reader's attention to a single 
straight-line-element, along which sources of constant strength 
are distributed. This simple case allows us to explain the 
basics of the panel technique. The source strength A is de- 
fined as the volume of fluid discharged per unit area. Since 
fluid is ejected perpendicular to the panel's surface in both 
directions, discharge velocities are half of the source 


strength. The boundary condition for an inclined panel requires 
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that the normal component of the free stream velocity is 


balanced by this discharge velocity (see Figure 3.1). 


== VY В (23224117) 
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Figure 3.1. Boundary Condition at One Inclined»Panel 


This relation between a geometric quantity (6) and the unknown 
source strength establishes tangential flow on the panel 
ъщигасе. 

Things, which are obvious for a single panel, become 
Slightly more complicated for a structure of panels. Mutual 
interference of source panels, i.e., each panel induces a 
velocity at other panels, causes the complication. While the 
boundary condition of a single panel had been set up by glancing 
at a simple geometry sketch, we now better switch to a syste- 
Matic procedure, emphasizing the concepts of velocity poten- 


tial and superposition. We consider a 2-D closed body, 


Jl 


approximated by several panels and inclined at an anale to 
the oncoming flow. Our goal is to derive a relation for the 
unknown source strengths from the condition of tangential flow 
at the control points. To this end we will give the velocity 
potentials of a single source, a source panel, and a closed 
body built up by a source panel, in that order. 

Radial streamlines and concentric circular equipotentials 
characterize one of the very basic potential flows, the single 


source flow. Its velocity potential is defined by 


Single source ф(х,у) = s+ ine (3.2) 
where r is the distance from (x,y) to the source. Arranging 


Single sources on a straight line element corresponds in 
terms of the velocity potential to a summation of single poten- 
tials, in the limit of an infinite number of sources to an 
integration over the panel length. Thus the velocity potential 
Of a source panel can be written as 

0, 


source panel $(x,y) = = Í Ш.Е "Сс Mc 
0 


m source panels, representing a body, induce a flow field, 


whose velocity potential at any point (x,y) is given by 


Ха m 
$ (x,y) o NE г. 48. (3.4) 


|| 
И 33 
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И mmg 


1: 128037) = 
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( хв; ‚Ув; ) Control points are 


midpoints of panels 






Boundary points 


Во. Designations for Calculation 


We cali $ the potential of the flow disturbances due to the dis- 


Beaeemene flow, The total velocity potential of a nonlifting 
flow results from a superposition of this displacement flow 


and a uniform flow, which is inclined at an angle a to the 


вахте. 
Т 21 
O(x,y) = У, созах + У 51п су + 1 27 | En Е ыг 


Recalling the definition of the velocity potential (velocity 


15 Gradient Of potential), the boundary condition of tan- 


gential flow takes the form 


Q2 


сша 0 on the surface 


Q? 


Applying this condition within the framework of the panel 


method provides a system of equations, 
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3 | ' 8 ( d 
т--(іп 1..)45. = 

2 i ) J 

fOr 20 (3.25) 


which establishes zero normal velocity at all control points. 
This linear system can be solved for the unknown source 
strengths after the integrals have been evaluated. 
Concept of influence coefficients 
The influence coefficient, СЕС is defined as the normal 
velocity at the 144 panel due to a source distribution of 


27-strength at the fe panel. 


Is, | TEE ri) ds, (3D 


The contribution to the normal velocity at the qum panel by 
the actual source distribution of thg Del panel is the product 


Ou 12/2 and the influence coefficient Iij- TO compute the 


influence coefficient we must substitute 


/ 2 2 
E = (x 2х. + Ум, =e 


1) М; 


in Eq. (3.6) and сакгу ойк тестте еш алш е 


Е хи-хи а 
- | E 11 D ишы (3.7) 
Фи ха 8 (Ум -У-) 
4: 1 
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where 





ах. 

дп, = Ш 

ду, 

SE = sin б. 
i а 

X. = x + S. cos 6. 
J В. J J 
Е = + S. Sin ð. 

5 ^B. j j 


The integration covers the length of the aoe panel. Еяна!1у 


we obtain 


- (x ven Coso. Ісоз8. + - +S.siné.) 5108, 
) pue Ug thu 


j М. Be 
i) F J ВИ ы AS кь _ E 
0 Dp a 2 ) ] ^ ey. 27227108 ) 1 


(2220) 


Equation (3.8) is expressed in terms of 5. only and, after 


some manipulations, the integral mav be written in the form 


1 b -cs. 
E aM алы (3.9) 
1J Ш 220011. 
J J 
where 
b = (Хи “Xp ) COS Bi * (Ум -Yg )Sin В. 
1 1 1 

с = cos шон В; + sin ln В; 


Je 


ИВ 


е Эзе -х ) iod 5119. (Ум -Ув ) 
1. 5 1 


Е 2 - 2 
Ё = (XQ 7Xg ) + (Ум Yp.) 
1 J 1 J 


The integration can now be performed analytically, 


1. 25.-е 1. 
.+f | ] о 2 — 


] 
j 
O утгаа 


© 2 
Е 5140 155-еѕ (3.109 


ЕЕ 
1) 


Determination of unknown source strengths à 
Equation (3.5), expressed in terms of Eq. (3.9) and divided 


DY Va, takes ti E orm 


TAL + 'Т.. = - cos acos B.- 51105118. (3.14 
i 1) T 1 


(1-4. 
M dg» 053 
H- 

> 


where 
J 
2пу 


or in the more convenient matrix ст 


! от: -D 1 1 

Iii Ij; ч К Т Ay COSQCOSP| sinasin, 
! чив = 1 1 

1-4 А2 созосозв, 51005116, 
! = =» 1 1 

i 1 Т А соѕасоѕВ -ѕіпаѕіпВв | 


(322228 
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The above set of linear equations can be solved for RS by 
Gauss' elimination method or any other linear equation 


algorithm. 


On-body velocities 


The velocity at the midpoint of the to 


panel, Ус сап 
qi 


be obtained by a spatial derivative of the velocity potential 


in tangential direction. 





Е 2 20 (х. ,у.) 
©. aS 
1 
ох ду m Е а, 
E 1 3E ! 44) 
Б Еее 1 хэ, "cs 45.) 
1 т 1-1 1 al 
(3-13) 
where 
OX. 
та = cos 0 
512 a 
1 
ду. 
Qc cc ү 
1 
Therefore 
Ус m 
Е = | | ‚+ оте Sie 4 
(0 ) cos acos 9, + sin asin m 4 | J ij ( ) 
сос 1 1-1 
where J.. = | Вт (аа цас denotes the tangential velocity 
1 1 
at the ы panel due to a source distribution of 27-strength 
at the Bn panel. 
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The calculation oT Ji; follows the same procedure as that 


ОЕ І о 
23 


5. b scs: 
з. =] ————L— às, (3-45) 
J Ооо ЕЕ J 
J J 
where 
b = (Хм -Хв ) Соз Ө. T (Ум "Урд ) ѕ1п 9, 
1 1 1 
с = cos б. cos 02251092 ти 
J 1 J 1 
е = 2 СЕТЕ + 2 sin. (Ym, Yp’ 
Qi J 1 
2 2 
Е ‘у usc ee 
1 2 1 J 


Positive signs of on-body velocities indicate that velocities 
are oriented in the direction of the surface coordinates, 
while negative signs indicate opposite directions of veloci- 
ties and surface coordinates. The positive direction of 
Surface coordinates is defined clockwise. Therefore positive 
values of on-body velocities are to be expected on the upper 
surface, negative values on the lower surface. 

Off-body velocities 

Streamlines can be determined by computing velocities at 


off-body points and using a numerical quadrature to progress 
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from one point to another point on a streamline. The velocity 


components can be expressed according to 


Сх. А) m à 
0 (х,,у:) a V,cos o* ) 27 Ира PE ds 
1 1-1 j i 
9%(х.,у.) 7 A: f 5 
DX v.) = = V sin а+ де ——- Lnr.. dS. 
ды =. ду. id 2T ; ду, 1-2 J 


Normalizing the above equations by the free stream velocity 


and abbreviating the integrals simplify the relations to 


У.) m 
= oso UT A! ES. (3.16) 
Ме. 1-1 35 
те 
ENEMY.) m 
= = sina+ ) А! т?. ТЕРТІМА 
Me. 3-1 шиг! 


x 


ШЕ апа 11. are again influence coefficients, whose evaluation 


can be adopted from the already introduced procedure. 


P. 
1 [pc e ETE S. 
5 as. (3.18) 
х 0 52 е5 +Е 2 
j j 
*j  b,-C,.S. 
тї. = f° ++ as, ЕРЕ) 
29) 0 Si -eS, +f : 
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where 


x I B: Y DE 
J J 
€ = | Em! = | 
Е cos us y sın 9s 
е = 2((х.-хь )соз zx + (Y;7yg )sin oe 
J J 
_ E 2 Е 2 
го е а Л На Yi 


C. CHERCULATORY ВВ 

While inviscid 2-D flow theories are not capable of pre- 
dicting drag characteristics, information about lift can be 
provided by them. Creation of lift is closely related to a 
type of flow called circulatory flow. We mentioned already 
that the flow around a lifting airfoil can be decomposed into 
two elementary flows, i.e., displacement flow and circulatory 
flow. Circulation and circulatory flows are the subjects of 
this section. 

The early approaches of airfoil theory emphasized a flow 
model in which the airfoil was represented by an infinitely 
thin vortex sheet only. This so-called thin airfoil theory 
predicts lift quite well, because lift depends primarily on 
circulatory flow. Unfortunately a straightforward extension 
of vortex sheets to "surface singularity" method is impossible. 


Therefore aerodynamicists have proposed a couple of flow models, 
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which allow the implementation of circulatory flows in 
"surface-singularity" methods. Examples are: 


(1) Smith and Hess represent circulatory flows by a com- 


bination of source and vortex distributions.  [Ref. 1] 
(2) Martensen prefers vortex distributions only, but states 
the problems in terms of the stream function. [Ref. 2] 
(3) Davenport makes use of linearly varying vortex 
distributions.» Ref. 3] 


Our approach follows the ideas of Smith and Hess. These circu- 
latory flows are composed of a vortex distribution, which is 
constant along all and for all panels, and a source distribution 
of conventional shape. 

We start at the very beginnings of vortex flows. Concentric 
circular streamlines and radial equipotentials characterize the 
flow field of a single vortex. Its velocity potential can be 


written as 


a20) 





L Y Yy 
single vortex ф(х,у) = ~- 5 7 arctan z3 


with (XY) as the center of the vortex. A structure of m 
vortex panels induces a flow field, whose velocity potential 


ШЕНІ роіпе (х,у) іс атуеп Ру 





m 1 1 Р 
ф (х,у) = Г 2. (- 217), arctan xx. аз. (2752189) 


ЕТ» field differs in two important points from the non- 
lifting flow field: 
(ПОМ с1о1а:є сс the condition of tangential flow. 


(2) The unknown singularity strength Il cannot be 
determined immediately. 
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The task of determining circulation must be postponed to the 
implementation of the Kutta condition. Temporarily we set 

the vortex strength equal to one. Tangential flow must be 
established by the aid of an additional source distribution. 
Strengths of this additional source distribution must be con- 
puted according to the condition that the normal velocities 
due to the vortex distributions at the control points are 
balanced by the normal velocities due to the additional source 


ДТС ПЛ 


Та а m — VN лг 
) = J д. шиг dS. = || шй : В 
ха цай т 1) 1 лае со пие SM " 
EIL шин: = Y 13 
Abbreviating the integrals by the above defined influence 


coefficients, we get 


га 
Wes E MENSEM" (3.25) 


1 J ij 


| wmm 3 


j 


where 


Nus are the unknown strengths whose effect is intended 
J to balance the normal velocities induced by a unit 
vortex distribution. 


га is the normal influence coefficient due to a source 
J distribution. 


TUR is the normal influence coefficient due tova vortex 
J distribution. 


x is the tangential influence coefficient due to a 
J source distribution. 


nd is the tangential influence coefficient due to a 
17 | 
уот тех eb M 
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Samce intluence coefficients of source and vortex distributions 
are related by е” = — the above equation can be expressed 
according to 


In 
Ш ГЭ - OS! for ix lm (5.28) 
= = 


É a ER 


By solving this system for ae we determine the properties 


emer culatory flow of unit strength. 


Calculation of disturbance velocity due to unit circulatory 
flow 

The disturbance velocity, и. is composed of two parts, 

Бие Ое со the constant vortex distribution, the other due to 


cene additional source distribution 


In In 
QUEE Nc E MESES (3.24) 
D 


Making use again of the relation between influence coefficients 


(119) = 1(3),, we have 
15 Yj 
m 
Идеи Са маи ии I (3.25) 
i ОВИЕ 3 5 


В. ЕЕ Т7:МС А COMBINED FLOW 


THe Kutta condition serves as matching condition for nonlift- 


ШИН па circulatory flow. These two basic flows must be 
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Superimposed such that flows of upper and lower surface merge 
smoothly at the trailing edge. This crigqimaly versione ote 7] 
Kutta condition is usually substituted by the condition of 
zero load (or equal velocities on both upper and lower sur- 


face) at the trailing edge (see Figure 3.3). 


Figure 3.3. Single Airfoil: Superposition of Nonliteiag 
and Circulatory Flow, Controlled by the 
Kutta Condition 


Since the panel method does not permit the evaluation of 
velocities at the trailing edge, the Kutta condition is satis- 
fied approximately by requiring that velocities at the control 
points of the rearmost panels have equal magnitude. Therefore 
the rearmost panels should be chosen short so that flow at 
their midpoints will effectively represent that at the trailing 
edge. 

Determination of circulation 


ЕЙ 


Suppose pee and m panels are the closest panels to the trail- 


ing edge on the lower and upper surface, then we can write the 
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Mircea condition as 


-у (м) м ы е v) + та (3.26) 
where 99) denotes the tangential velocity in nonlifting flow. 
Нешев топ (3.26) can be Solved for the’ ctrctlation T. 
Calculation of on-body velocities and of pressure coefficients 
Three parts contribute to the total velocity: free stream, 
disturbance due to displacement flow and disturbance due to 


ligo flow. бау 2. designates the velocity due to the 


nonlifting flow including the free stream component and vv) 
represents the velocity due to a lifting flow of unit circu- 


lation. Then the total tangential velocity at the midpoint of 


the P” panel is given by 


Еа) амал UP Caer 
1 ШІ al 

Once the velocity has been computed, the pressure, customarily 
expressed by means of a dimensionless coefficient Б 15 


determined by Bernoulli's equation: 


Pi Po Me? 
ay ИЖ ^5 5 los (с) (57720) 
L 5 РУ. 85 


Addendum: Моге than one body configuration. 
One of the main advantages of the panel method is its 


easy extension to multi-element airfoils. As a matter of 


45 


fact even flow past an infinite number of bodies can be solved 
by means of the panel method, if these bodies are arranged in 
form of cascades. The minor changes, which are necessary to 
apply the panel method to a finite number Of bodies eine єє. 

(1) The overall scheme must provide a circulatory flow 
for each lifting body. (The number of nonlifting 
flows remains one.) 

(2) Flow past each body with lift is subject to a Kutta 
condition. Accordingly the numbers of equations 
requiring zero load at the trailing edge equals the 
number of circulatory flows, which allows the 


definite determination of each lifting body's 
circulation: 


Figure 3.4 illustrates the superposition of nonlifting. and 


circulatory flows for a two element airfoil. 


Figure 3.4. Two Element Airfoil: Superposition of 
Nonlifting and Circulatory Flows 
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Е. EXAMPLES 
This section illustrates the capabilities of the program 
"PANEL" which can be applied to 2-D potential flow problems 
past one or two bodies. 
miewepast Ore™=circular cylinder 
Шис сос panci technique is applied towthe flow past 
a circular cylinder. This case is regarded as nonlifting, 
i.e., the cylinder experiences no force perpendicular to the free 
ЕЕ.  ASsOsketched іп Figure 3.5, the surface of the cylinder 
is approximated by eight panels of equal width. For zero angle 


pL tack, Eg. (5.5) reduces to 


À 


i DUE 5-2.) 
v cos; + 5 + 


À 
— ee. с. 
2 


т. s 0 (22129) 
1 1 Í 
i 


1 сезі 


3 
J 


Solving a set of 8 simultaneous algebraic equations, the source 


strengths and the pressure coefficients can be determined. 





Це 
ee 
Figure 3.5. Arrangement of Panels on a Circular Cylinder 
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The results are shown in Figure 3.6 where they are com- 


pared with analytical results (Ср = 1 СУ 


Феееецеевееееееееееееоееее еееейееоееееоео бефоееееееоеееФтееееееееееейдееееееееегеее 
9 


ФеесееоеефеееоеееееФееоееееейеФәеКеЭееегооеоОФееееефооееооооооеФФеезейоеееоооефооеееге 


СР 
-1.0 





0.0 90.0 - 180.0 2700 360.0 
ТНЕТА 


Figure 3.6. Pressure Coefficient on a Circular Cylinder 
Obtained by Using Eight Source Panels 
(Marked by 0) in Comparison with the Exact 
Ec OI 


This example demonstrates the power of the panel method. 
However the reader should be aware that only 8 panels are not 
sufficient to describe the geometry in most of the cases. 
Basically the achieved accuracy depends on both the shape of 
the body and the panel configuration (number of panels and 
local widths). A closer spacing is advisable in regions where 
severe changes of the pressure distribution are expected 


(e.g., leading edge). 
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еи расе a pair of circular cylinders 

Two circular cylinders are arrangedside-by-side in a uni- 
form stream. The surface of each cylinder is replaced by 50 
panels of equal width (see Figure 3.7(a)). The computed 
p eene distribution on one of the cylinders is shown in 
ВЕЧЕ 3.36). The reader shall pay some attention to a 
comparison between the flow past one and the flow past a pair 
of cylinders. Obviously the maximum velocity is increased 
by the existence of a second cylinder. The closer the two 
cylinders are arranged, the higher the maximum velocity. 
While the stagnation points in a single cylinder flow are 
located at the farthest down and upstream points of the cylinder, 
the disturbance of a second cylinder causes the stagnation 
points to move towards the other cylinder. The streamline 
рше Еше, given in Figure 3.8, should provide a deeper under- 
ана ОЕ this kind of flow. 

Hie. past two element airfoil 

The main goal of leading and trailing edge devices is to 
obtain a higher lift coefficient. We will investigate the 
effect of a single slotted flap on the pressure distribution 
of the main airfoil. 

[fe pressure distributions of a single airfoil and of an 
ЕШ ОК Етар Combination are Compared in Figure 3.9. Тһе 
coordinates of both main airfoil (a NACA 4412) and flap are 
listed in Section F (sample input data). The results indicate 
that lift increases more than 50% by using a slotted, 21.5 


degree deflected flap. 
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50 panels 


-90.0 -60.0 -30.0 0.0 30.0 60.0 90.0 


ТНЕТА 
(b) 
Figure 3.7. (a)Arrangement of Panels on Two Cylinders 


Side by Side 

(b)Calculated Velocity Distribution on One 
of Two Identical Circular Gy landers Whose 
Centers Are One andma Оше a amemems 
Apart 
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~20.0 ~ 25.6, 


- 15.0 
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Figure 3.9. Pressure Distributions on a Single and 
a Flapped Airfoil 


Аіг ҒОЛ о ле Мете 
Flow past an airfoil in ground effect is another ар рї са 
tion ОГ єг ээ 20 
The boundary condition at the ground requires vanishing 
normal velocity there. We meet indirectly this condition or 
arranging the second, imaginary airfoil such that the ground 


becomes an axis of symmetry of this "new" flow field (see 


De 


Figure 3.10). Since an axis of symmetry must be impermeable 

to fluid particles, the desired flow is obtained without 
emplicitly satisfying the boundary condition at the ground. 
This kind of flow is a challenge to aerodynamicists for several 
reasons. Whenever an airplane takes off and lands, it passes 

a zone where flow is severely affected by the proximity of the 
ground. Wind tunnel experiments must be corrected for wall- 
effects, quite a similar situation with grounds below and above 
the airfoil. And there was a German experimental seaplane 

that makes use of flying very close to the sea level. However, 
our numerical experiments will tell only one part of the story 


because all these flows are highly 3-dimensional. 


Peal 


<< 22--20- 
м. 
Ground . 


GU Quo 
..... 
а 
ооо 
“ае а 
4710026 S 


ЕЕЕ 102 АТ ий іп Ground Effect 


Let's first question how does the pressure distribution 


СИЛ Ее спе ground. Figure 3.11 shows that lift is reduced 


53 





80 01 0% 03 04 05 06 07 08 00 19 


Figure 3.11. Pressure Distributions on a Single NACA 4412 
and on a NACA 4412 in Ground Effect 
(ИЕ ЕО. 2. СОУ) 
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Figure 3.12. The Lift of a NACA 4412 Near the Ground 
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оп the upper surface and increased on the lower surface. Іп 
the particular case the overall lift gain is about 15% of the 
Шин tire cree air, bul we Might not always expect a lift gain. 
The actual balance between lift reduction on the upper and 
lift increase on the lower surface depends on both distance 
from the ground and angle of incidence. Figure 3.12 confirms 
that there are cases with less lift than in free air. High 
angles of attack and moderate distances from ground are sus- 


ceptible constellations to lift loss. 


Е. I/O--DESCRIPTION AND LISTING OF THE PROGRAM "PANEL" 
This program calculates non-lifting and lifting potential 
flow past one or two bodies. Any 2-dimensional shape and 
any angle of attack, which do not cause flow separation, 
are acceptable. 
Input data 
The data must be arranged in the following order: 
(Т) Header card; 


(2) Coordinates of first body cards (variable number of 
cards); 


Second body control card; 


(4) Coordinates of second body cards (variable number of 
Cards) 


Items 3 and 4 are used only for the 2-body case. The actual 
instructions are as follows. 
Header card 
l-10: Number of bodies (integer) 
11-20: Number of points of the first body (integer) 
21-30: Angle of attack in degrees (real) 
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Coordinates of first body cards 
The input procedure of body coordinates requires the follow- 

ing sequence: start at the trailing edge, progress on the 
lower surface to the leading edge, return on the upper surface 
to the trailing edge and finish with the trailing edge. The 
trailing edge of closed bodies is input twice, as first and 
last point. However airfoils with finite trailing edge thick- 
ness should be treated as non-closed bodies, i.e., the last 
point input is not the first point repeated. 

1-10: (X coordinates of the points defining the body (real) 
11-20: Y coordinates of the points defining the body (real) 
Second body control card 

1-10: Number of points of the second body 
Coordinates of second body card 


The X- and Y-coordinates of the nd 


body are input in the 
same format as the coordinates of the first body. 
Output 
There are two kinds of solutions, non-lifting and lifting, 
both of which are preceded by the following column header. 
PANEL Х Ү V Є 


p 
where 


PANEL is the number of the panel; 


X and Y are the coordinates of control points (not 
boundary points); 


V denotes the relative тео он ое 


C denotes the pressure coefficient. 
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Sample problem 
This sample illustrates program input and output. The 
data refer to the airfoil-flap example of Section III.E 


(see Figure 3.9). 


Ээ! 


Input 


н н нә а на rm 
е е Ф е е е е е ж ө 


ктк о оа о ыт коа USD UTI 
© 
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N 
о 


каак таа т | 
e 
Un 
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Output 
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онн ОО О О М М № № № № СЧС № + © © 


52361 
47182 


55519 
22992 
07829 


ЭО ООЭ ОО ОО ООС«ОООСОО 


, 
кз 
o 


оу ой 
ООО ООС«ОО оным лм) (а ч са Ал СХ Со го 


(11 
оон № 


СР 


219955 
. 34430 
.49083 
.566 38 
.60157 
.6 34532 
-67897 
272585 
77759 
.81548 
287557 
.94714 
499587 
.99041 
«77108 
.19848 
.89880 
.97764 
"54344 
‚37716 
„15257 
‚50686 
„ВИЗИ: 
„3536759 
‚02855 
. 78674 
. 30540 
‚87199 
‚56299 
. 29451 
. 05920 
.78706 
.66265 
.19956 
.25482 
257509 
.81064 
.85010 
.68455 
.24862 
252159 
.25482 


Program listing 


 СЕ-Есесгссссссввссесесвсгесссссссссссссссссссссссссссссссссссссссссссс 


THIS PROGRAM CALCULATES 2-D POTENTIAL FLOW PAST 1 OR 2-BODIES 
AT ANY ANGLE OF ATTACK. 


WRITTEN BY : CAPT LEE,HEE WOO 
DATE : NOV.28 1985 


NOTE : MAXIMUM NUMBER OF PANELS = 200 


БЕРСЕ ЕРЕССЗЕРЕСЕРЕССССССССЕОРОЕСОССГСССССССССССССССССССССССССССССС 
ОТМЕНЪТОМ 2200,200), ХВС200), УВС200), ВЕС200),ТНС200), МС200 ) 
DIMENSION WKAREAC65000),XMC200),YMC200),VVVC200),VCC200) 
DIMENSION CPC2000,VVC20050,YC200,2002,AC2005,V1C2002,31C200) 
DIMENSION VC1C2000,VC2(2000,V2C2000, VTC200) 


C 
Ca READ INPUT DATA CFOR FIRST BODY) --- 
Е 


КЕАБ( 4,1) МВ, ММ, АМ 
1 FORMATC2I10,F10.5) 
АМ = АНХ5.141592/180. 
0 10 Т = 1, ММ 
READC4,110 XBCIO,YBCID 
11 FORMAT(2F10.5) 
10 CONTINUE 


9 ООО ООООСОО 
ООООООО 


С 
C --- CALCULATE MID-POINIS OF PANELS AND ANGLES THETA --- 
С 
М = ММ-1 
DO 12 I=1,N 
= 1+1 
ХМ(1) = (ХВ(1)+ХВ(К))^/2. 
ACID = (YBL CIDFYBCK))/2 


THH = (YBCK)-YBCI))/C(XBCK)-XBCI)) 

THCI)= ATANCTHH) 

IFCXBCK).LT.XBCI)) THCIJ=THCI)+3.141592 
12 CONTINUE 


C 
EL = CALCULATE PANEL LENGTHS --=- 
C 
DO 13 I = 1,N 
K = I+l 


ACI) = SQRTCOXBCK)-XBCI) )¥¥2+0CYBCKI-YBCI) ) ¥¥2) 
13 CONTINUE 
IFCNB.NE.1) GO TO 600 
ESGGSSSSSGSOCOECCCEGBGCECCCCCCCCCCCCCCCEGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C C 


С NON-LIFTING PART C1-BODY) С 


© C 
Е ССС ЕСССССЕСССССССССССССессссссссссссссссеесессссссссссссссс 
DO 14 1 = 1,N 


C 
С --- CALCULATE ANGLE BETA AND NORMAL COMPONENT OF FREE STREAM 
С VELOCITY ==- 
C 
BEGIO- HHCTI+5. 14159272. 
УСТ» = -CCOSCBECI)DIXCOSCANIF+SINCBECIDIXSINCAND) D 
DOTIS J = 1,N 
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=== CALCULATE INFLUENCE COEFFICIENTS OF NORMAL VELOCITY === 


IF Cl. EQ) СО ie 
(ХМС17:-38()))ХхС05ОВЕСТЭЭ «СҮМСТЭ-Ү32 ЕЕ ВЕС 
COSCTHCJ) )XCOSCBECI)DI+SINCTHC J) XSINCBECI)) 

COSCTHCJOOXCXMCIO-XBCJOO*2.XSINCTHCJOOXxCYMCIO-YBCJO) 
CXMCIO-XBCJOOXX2*CYMCIDO-YBCJOOxx2 
SQRTC4.XF-EXX2) 

2(1,4) = TEGCB,C,E,F,H,ACJ)) 

15 CONTINUE 


Z(I,I)=3.141592 
14 CONTINUE 


OO’ 


ии ж ии 


INNO 


С | 
С === SOLVE SET OF LINEAR EQUATIONS FOR SOURCE STRENGTHS --- 
E 
IDGT = 0 
CALL LEQT2F (Z,1,N,200,V, IDGT, HKAREA, IER) 
C 
C --- CALCULATE INFLUENCE COEFFICIENTS OF TANGENTIAE УЕ 0 08 2--- 
С 
DO 16 I = 1,N 
DO 17 J = 1,N 
ТЕСТ.Е@ ШУ СО ТО 1? 
В = (XMCID)-XBCJOOXCOSCTHCIOO*CYMCIO-YBCJOOXSINCTHCID) 
С - COSCTHCJODOXxCOSCTHCIOOTSINGDHCJOOXSTINCDHQIMED 
Е = 2.XCOSCTHCJO)OXCXMCIO-XBCJO2) *2.XSINCTHCJODOXCYMCIO-YBCJ2) 
Е = (XMCID-XBCJ2)2XX2* CYMCID) -YBCJOOXx2 
H * SQRTC6.XF-EXX2) 
YCI,J) * TEGCB,C, E, FH, ACJ2) 
17 CONTINUE 
ҮСІ,1)=0.0 
16 CONTINUE 
С 
С --- CALCULATE TOTAL VELOCITY AND СР AT MIDPOINTS OF EACH PANEL --- 
С 
WRITEC8,95) 


95  FORMATC//7,25X, 'NONLIFTING SOLUTION',/77 
5105. PANET АТАНЫ ло ЭЛ 
DO 18 I = 1,N 


5=0. 
DO 19 J = 1,N 
S = S+V(J)XY(I,J) 
19 CONTINUE 
УМСІ) з СОЗСТНС125 )ХСО5САМ) 451МСТН(С1) )Х51ЇМСАН) 455 
ЕРСІ) - 15-ууфстрожка 


WRITEC8,95) I,XMCIO,YMCIO,VVCID,CPCID 
95 ҒОКМАТ(10Х,15,5Х,2Ғ10.5,5Х,2Ғ10.5) 
18 CONTINUE 
сссссссссссссссессссссссессессссесеЕсеЕссеЕсеРЕЕЕЕРЕ сг ЕЕЕ ЕЕ МЕРССССЕСЕИЕ 
> С 


C LIFTING PART C1-BODY) C 
> С 
CCCCCCCCCCCCCECCCCCCCCCCCCCCCCCCECECCOECPERECCCCCEESCOCCGECCDOEPECCCCCEOIP OMS 
[E 


C 
С --- CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 
C STRENGTH === 
C 
DO 20 I - L,N 
5 = 0. 
DO 21 J = L,N 
= 3tYCI,J) 


él CONTINUE 


20 CONTINUE 
IDGT = 0 
CALL LEQT2F (Z,1,N,200,V1, IDGT,WKAREA, IER) 


VICID=S 
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ооо О 


ООО ООО 


С 


96 


24 
25 


C 
600 


GOOG 


ооо 


31 


32 


23 


34 


35 


CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOW OF UNIT 


STRENGTH === 
DO 22 I = 1,N 
S = 0. 
DO 23 J = 1,N 
S = S*YCI,JOXV1CJO*ZCI,J) 
CONTINUE 
УССТ) =5 
CONTINUE 


CALCULATE VORTEX STRENGTHS BY KUTTA CONDITION --- 
SI * -CVVCIO*VVCNDD/CVCCIO*VCCNDD 
CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL --- 


WRITEC8,96) 
FORMAT(C///,27X, 'LIFTING SOLUTION'!,// 


X10X, 'PANEL',5X, 'XM' ,8X, 'YM' , 11X, 'V' , 10X, 'CP') 


DO 24 I =1,N 
VVVCID = УУ(Т) + 5ТХУС(Т) 
CPCI) = 1.-VVV(CIOxx?2 
WRITEC8,93)1, XMCIO,YMCIO,VVVCIDO,CPCID 
CONTINUE 
мЕОКНАТСЭРИ ЭЭ 
GO TO 700 


READ INPUT DATA CFOR SECOND BODY) --- 


READC4,31) MM 
FORMATCI10) 
DO 32 I = NN+1,NN+MM 
READCG,11)XBCI),YBCI) 
CONTINUE 
М = NN*MM-2 


CALCULATE MID-POINTS OF PANELS AND ANGLES THETA --- 


DO 33 I=N+1,M 

К = 1+1 

ХМС1) < (ХВ(К) +ХВ(К+1))/2. 

XMOI) с СҮБСКЭТҮВСКТІ))/2. 

THH = CYBCK+1)-YBCK)I/CXBCK4+1)-XBCK) ) 

ТНСГ) = АТАНСТНН) 

IFCXBCK+1).LT.XBCK)) ТН(Т) -ТН(Т) +5.191592 
CONTINUE 


CALCULATE PANEL LENGTHS --- 
DO 34 I = №+1,М 


K = [+1 
ACI) = SQRT((XBCK+1)-XBC(K))XX2+(YB(K+1)-YB(K))XX2) 
CONTINUE 
DO 35 I = М+1,М+1 
XBCI) = XBCI+1) 
ЕСІ) = ҮВ(Ї+1) 
CONTINUE 
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ОО С) 


O0 0€0 


ОСО 


ОО С 


OOO 


оо 


37 
36 


39 
38 


41 


40 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCE CE CECE CCC CCC CCE CCE CCCCCEOCCCECECC CEC CCCCeeee 
С 


NON-LIFTING PART (2-ВООТЕ5 ) С 


С 
CCCCCCCCCCCCCCCCCCCCCCCE ECC CC CCC CC CEC EC EE CCC CC CCC CCC CCCCEC CCE Gere CCC eer 


DO 36 I = 1,M 


CALCULATE ANGLE BETA AND NORMAL COMPONENT OF FREE STREAM 
МЕРОСИ s: 


ВЕСІ)- ТІПСІ)Ғ5 141592727 
VCI) * -(COSCBECIOOXCOSCAND*SINCBECIOOXSINCAN)) 


CALCULATE INFLUENCE COEFFICIENTS OF NORMAL VELOCITY --- 


DO 37 J = L,M 
IFCI EQ 202 CO 1017 


В = (XMCIO-XBCJOOXCOSCBECIOO-*CYMCIO-YBCJOOXSINCBECID) 
C = COSC THCJ))XCOSC BECI))+SINCTHC J) DXSINC BECT)DD 
Е = 2.XCOSCTHCJO99XCXMCIDO-XBCJ29*2.XSINCTHCJOOXxXCYMCIDO-YBCJ2) 
Е = (ХМСІ)-ХВеЈ) ) жж2+(ҮМ(І) -ҮВ(Ј) )хх2 
Н з 50КТСФ4.ХР-Ехх2) 
ZCI J) =T EGB. C.E, F; H ACIDI 
CONTINUE 
2(І,І2-5.141592 
CONTINUE 


SOLVE SET OF LINEAR EQUATIONS FOR SOURCE STRENGTHS --- 


ТОСТ = 0. 
CALL LEQT2F (Z,1,M,200,V, IDGT, HKAREA, IER) 


CALCULATE INFLUENCE COEFFICIENTS OF TANGENTIAL VELOCITY --- 
DO 38 I = 1,M 


DO 39 J = 1,M 
IFC(I.EQ. J>) GO TO 39 
B = (XMCIO-XBCJOOXCOSCTHCIOO*CYMCIO-YBCJOOXSINCTHCIO) 
C = СО5(0ТН(С.) )хСО5(ТН(1))451МСТНС.,Э ЭХ51МСТН(1)) 
Е = 2.X€05C THC) >*XCXMCID-XBCJ) Э-27351М 11619 уУх СО Т»-ҮВС/8) 
Р = (ХМ(С19)-ХВС/2 Эхх24(ҮМ(1)-ҮВС)))хх2 
H - SQRT(C4.XF-EXX2) 
YCI,J) * TEGCB, C, E, FH, ACJOD 
CONTINUE 
YCI,1I)=0.0 
CONTINUE 
WRITEC8,95) 


CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL --- 
DO 40 I = 1,M 


5=0. 

DO 41 J = 1,М 
9 - StVCJOXYCI,J) 

CONTINUE 
VVCID) 7 COSCTHCIOOXCOSCANO*SINCTHCIODXSINCANO«*S 
CPCI) = I.-VVCTIOXX2 

WRITEC8, 93) I,XMCID),YMCI)J,VVCI), CPCI) 
CONTINUE 
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ы с = сс ссссссессссссеБЄСесССссСссС 
С 


С 


LIFTING PART (2-BODIES) С 


С С 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


ОООО 


ООО 


ООО 


ОООО 


43 
G2 


45 
44 


47 
G6 


49 
48 


CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 
STRENGTHAFASTTFIRST BODY —-- 


DO 42 I = 1,М 
S= 0. 
DO 43 J = 1,М 
G - 1l. 
тв. СТ.М) 050. 
$ = StYCI,J)%G 
CONTINUE 
У1(1)55 
CONTINUE 
ТОСТ =0. 


CALL LEQT2F (2,1,M,200,V1,IDGT,WKAREA, IER) 


CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOW OF UNIT 
STRENGTH PASTEUR S! BOD T=- 


DO 44 I = 1,M 
5 = 0. 
DO 45 J = L,M 
G = 1. 
PECI TGT. ND G30: 
S = S+Y(I,J)XVIC(J)+Z(I,J)XG 
CONTINUE 
ме S 
CONTINUE 


CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 
STRENGTH PAST SECOND BODY --- 


DO 46 I = 1,M 
э = 0; 
DO 47 J = 1,M 
G= T: 
ТЕС Дд БЕ МО ИВА 
S = ЭТҮСІ;22Х6 
CONTINUE 
V2CI)=5 
CONTINUE 
IDGT=0. 


CALL LEQT2F (Z,1,M,200,V2, IDGT, HKAREA, IER) 


CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOM OF UNIT 
STRENGTH PAST SECOND BODY --- 


DO 48 I = 1,M 
ous D. 
DO 49 J = L,M 
бо 1. 
ТР. ТЕ. Ц) 6-0. 
5 = S*YCI,JOXV2CJO*ZCI, Jo xG 
CONTINUE 
VC2(I)-5 
CONTINUE 
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ОО О 


50 
700 


ООО 


--- CALCULATE VORTEX STRENGTHS BY KUTTA CONDITION --- 


Х1 5 VC1C1)+VC1(0N) 
X2 = VC2C1)+VC20N) 
X3 = -VVC1I-VVOND 

XG = УС1(М+1) +УС1(М) 
X5 = УС2(М+1) +УСЕ(М) 
X6 5 -УУ(М417-УУ(М) 


512 з (Х3ХХ98-Х1ХХ6)/(Х9ХХ2-Х1ХХ5) 
- (Х5-Х2Ж5122/Х1 


CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL --- 


WRITEC8, 96) 
DO 50 I =1,M 
VTCI)D 5 УС1(01)х51144УС2(1)х512 
VVVCID = УУ(Т) +УТ(Т) 
CPCI) -S 1 = 
МКТТЕСЗ,93)1,ХМС1),ҮМСТ) ,УУУСТУ,СРСТ) 
CONTINUE 
WRITEC6 ,97) 
FORMATC1X, ‘COMPUTATION COMPLETED') 
STOP 
END 


--- THIS FUNTION EVALUATES THE INTEGRALS CINFLUENCE COEFFICIENTS Sa 


FUNCTION TEG(CB,C,P,Q,R, 3) 
ТЕКМ1 5 АШО0б((5хх2-Рх540)/0) 
ТЕКМ2 5 АТАМСС2.Х5-Р)/К)-АТАМС-Р/ЛК) 
ТЕО ә "СХТЕКМ1/2,%(2.хХ8-СХРОХТЕКМЕЕ 
КЕТУКМ 
END 
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ТУ. BOA METHOD 


Aw 6LILN@RODUCTION 

The thin shear layer equations are more complicated than 
Laplace's equation because they are nonlinear. This chapter 
presents the box-method, which can be applied to the solution 
of the thin shear layer equations. The box method was intro- 
duced by Keller in 1970 [Ref. 4]. 

One of the basic ideas of the box method is to write the 
governing system of equations in the form of a first-order 
system. This system is solved by finite-difference approxima- 
tions and Newton's method is applied to solve the equations. 
Finally, the resulting linear system is solved by the block- 


elimination method. 


В. FALKNER-SKAN TRANSFORMATION 
The thin shear layer equations for incompressible laminar 


flow take the form 


Qu oV 
мин — эь» d 
эх E ду 0 1) 
du du t dP 3^u 
sax C By 5 трах? у 2 иг 
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and at the edge of the boundary layer 


ТЕ 15 convenient to reformulate the equations using the 
streamfunction and the similarity concept. Therefore the 


Еа1Кпег-5Кап transformation is introduced. 





U 1722 ве: /? 
n ren! y = эн (4.4) 
(х,у) = (U, vx) I ^ (x , n) (4-45) 


Substituting Eqs. (4.4) and (4.5) into Eq. (4.2), we get the 


transformed momentum equation for 2-D laminar flows. 


Dn Ie е 242 E сох E 
Е + — ЕЕ [ЕЕ иг s = f x! (4.6) 
where 
dU 
m = x е 
О ах 


(dimensionless pressure-gradient) 


with the boundary conditions 


(4.7) 
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Шина Sea. function ofen only the right-hand terms of Eq. 
(Geom be zero. Then this will be a third-order ordinary 
differential equation whose solution is well-known as a 
Монг Elow. Būt, if f is a function of n and x (non= 


similar flows), we need a numerical procedure, such as the 


DOs method. 


See ERICAL FORMULATION (BOX METHOD) 
First of all, the coordinates (x,y) of a given geometry 


ese oe transformed to coordinates (¿,n) to apply the Box 


method (see Figure 4.1). 


> эсабпатіоп Potnt 


+» 
хас 


АТГ 011 





шинэ 
ТЭГ» 


НЫ... 








с 2 2 
51-1 э ЕО 


Transformed Ccordinates of Upper Surface 
Airfoil and Net Rectangle for Difference 


Approximations 


Figure 4.1. 
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The boundary layer thickness in transformed coordinates is 
nearly independent of the streamwise distance and can be 
represented by a fixed number of profile points at fixed 
spacing. 

One of the basic ideas of the Box method is to write the 


governing system of equations in the form of a first order 


system. We write Eq. (4.6) in terms of a first-order system 
of PDE's 

Е - а(Е6,п) (4.8а) 

u' = у(ё,П) (4.8b) 

(Ку)! + (EES) ev + КТЕ | = &(а Е - у 52) (4.8с) 


where a prime denotes differentiation with respect to п 


and ВЕЕ A 


with the boundary condition 


£(£,0) - 0,  u(£,0) 0, от (25 


We denote the net points shown in Figure 4.1 as 


And we can introduce the following approximations: 


A) Coordinates of mrzdpoints (t£ 177 1) and net functions 


1-5 1-5 
йг Stands for f, u or v) 2 2 
2 JT 
i-l 
2 х на m Se ee E 
EE = Е +9 ) E = 2(9- ШЕРТ (225012) 
2 


where [ І; елшеш Ее quantities (f or u or v) at point шиг 


В) Finite-difference approximation 
From Eqs. (4.8a) and (4.8b), the centered-difference 


derivatives are 


1 1. 
25723-1 і 
- шиг (4.11а) 
J су 
u} -u | 
j 1 -l _ i и (4.11b) 
j TT 
After introducing these approximations into Eq. (4.8) and 


rearranging (the known quantities are moved to the right hand 
side), we get the equation (4.12c) which is centered about the 
Болт (С 177 1) · This represents the relationship of quan- 


Mm 2 = = 


tities beteen еһе points of the box. 


О (22124) 
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1 _ _ 2 
ЕЕ Ul 1 5-(У D E 0 (4,125) 
ioe aL i 2”? 
Usar b SUE ]l/h, + ау (Еу) THE alu) 1 
272 J 
| J i-l_i 2 ша! 
- т 21 - P M EL = E (421207 
372 372 ог о ШЕРІ 
мһеге 
5 
1-2 m.+l 
"S ae ау = 5 +а, а = m. +а 
J 
| -1 
p o NE 2 E 1 
в ES (-L 1 +a[(fv) 1764 ) 12: -m 
D 175 2857 Jaz 
2 lo „у. 1-1 
ll MN 55 иж! +] E 
L 1 { А E (Е ) mii (u) 11} 
Jum J Да CS 


The last of the above equations is non-linear. Therefore we 


introduce Newton's method to solve this system. We set 


сэм” Р 8.7) " 6g 57) ЭЭ (4.13) 


where the superscript in parentheses is the iteration counter 


with initial condition 


ша 


(0) Е (0) те 

NO -0 i-1 (0) i-1 (0) i=] 
Е ЕЕЕ . * = 8 а === 

j Е. у E V Е (4.14) 
(0) _  ,i-1 (0) Оо 
- = у us =z | V5 VI 


EuEseSducing Eq. (4.13) to Eq. (4.12) and dropping the quadratic 


terms in gi, we get (superscripts i and n are dropped for 


simplicity) 
I 
2 В. 
Bre оо! до. + S DE ub Eos) 
5. 
оц. = ое = ко + ESSE = E а. ТБ) 
ae + DOE CEDE + т 
+ (5 ба. + (S Tu: = Ла qc 
(55) 69. + (56). 69. 1 (r3); ( 
where all terms are explained in Reference 5, with the 
peuncary Conditions 
Of = 0 бу = 0 бу = 0 (4.16) 


БР ЕСЕОСЕ ELIMINATION METHOD 
This is a very effective way to solve linear difference 


equations, discussed by Keller in 1974 [Ref. 6]. We write 


25 


Eq. (4.15): in a matrix-vector form 


Aé = г (441178 
where 
а 50 E 
П б] 2 
A = ` _ ` S | 6 = r = я 
ТІ ae 
E 50] rJ 
where 
TE) 0 1 1/2 0 
А, = 0 1 0 A, = (S3) 5 (55) - (51) - 
0 -1 “h,/2 0 1 0 
k L | 
1 Е 0 0 0 0 
З= Sa) с шар С = 0 0 0 < < J-l 
j (S3); (S5); (51). 5 l ази 
| 0 -1 prom /2 0 1 EE 
-1 шэг 0 
Ba = Su Se). 5). ju S 
j E Е ( 2/5 < 7 
0 0 0 
where 
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ЕТ 
J 


= = би 02:21 111 
Ove 
k Jj 
цас 0 
E: - ae TENET = Dom 0 
E (05) 0 
| L 
(ri); 
гу = (r4) 4 


According to Keller's block elimination method, we have to 


еее the matrix A. 


А = PxỌ (4.18) 


where 
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I is the identity matrix 


1 0 0 
т Фе 0 1 0 
0 0 1 


From Eq. (4.18), we find that 


90 - А, (4.19а) 

Р.О. = В. 15212225 ? 4,196 
393-1 3 j J ( ) 
= A.-P.C. = 1,2, "ul 4.19c 

9. j j j-1 J ( ) 


Keller showed that the matrix 2 has the same structure as 


the matrix ЫГ From Eq. (4.19), we derive the elements of 
Р. 202012 
j an Q4 
(P31?3 РЛ ЕЕ 
Е = (52177 (22275 (Роз). 
0 0 0 
5 2 
(omnis штэ шээг 
0 zi REP 
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Pace lement оғ Р. апа 2 is explained in Reference 5. If 


we let 

oo + М (4.20) 
ШИРЕСЕЛІСІТІП Ее. (4.20) апа (4.18) іпео Еа. (4.17), 

РН = Е (4.21) 


шити Crom Ес. (4.21), we find that 


The elements of М. are listed in Reference 5. Finally, we 


deene Matrix form to Get d from Eq. (4.20). 


09 Со 90 Иб 
"NN 51 и, 
р Ex ан = (4.22) 
24122111 
Ср || ч "HT 


From Eq. (4.22), we find that 


0194 т Ws (4. 2за) 


127 


Du 44-16 
е J 


ра а 4.235 
394 ; 3 ( ) 


gni 


Therefore J can be obtained by calculating the terms 0., С. 


апа = (see Reference 5). 
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У. INTERACTION METHOD 


Аи» INTRODUCTION 

Interactive methods provide a special coupling between 
viscous and inviscid flows. They are intended to compute 
flows including separation. Thus these methods may be regarded 
as an alternative to the Navier-Stokes solvers, which are 
hardly engineering tools because of their huge computer time 
and storage requirements. 

The classical method to compute viscous flows past airfoils 
proceeds as follows: The velocity distribution is computed 
by any appropriate inviscid flow solver. The output of the 
inviscid flow solver becomes the input of the viscous flow 
о. ocOlvying for viscous flow Consists of integrating the 
boundary layer equations. Provided that the flow remains 
attached this procedure allows a reliable prediction of lift 
and drag. Information is transferred only once from inviscid 
to viscous regions. However, many flows cannot be modelled by 
one-time information transfers. | 

Separation bubbles and separated flows especially require 
a close coupling between viscous and inviscid regions.  Inter- 
action schemes provide a better exchange of information between 
these two regions. 

The elements of interaction schemes are: direct or inverse 


inviscid flow solver and direct or inverse viscous flow solver. 
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“Воцпаату Condition | 
Flow Direct Inverse 


Шах 1 (е 1101 * Zero normal * Prescription of 
velocity at velocity 
- the surface distribution 


Viscous - No slip . Мо slip condition 


condition 


e Prescription of 

- Prescription of displacement 
external thickness 
velocity 





The direct boundary layer method has the disadvantage that 
the boundary layer equations are singular at the point of 
separation. However, if the external velocity 15 computed 
by prescribing a displacement thickness (known as the inverse 
boundary layer method), they can be integrated through the 
point of separation. 

The next problem associated with the regions of reversed 
flow is numerical instability, because downstream marching 
procedures cannot be applied in regions of reversed flow. 

The most common approximation to get this instability under 
control, the so-called FLARE approximation, neglects the 
momentum transport term u du/dx in regions of reversed flow as 
long as this region is small. However, as the size of this 
region increases, the FLARE approximation becomes inaccurate. 
One of the procedures which can be taken into account is 
called the DUIT (Downstream, Upstream >Шсегастоа). 2221-1058 
of a sequence of alternating up and downstream sweeps. 

There are several types of recently developed interaction 


models. All procedures have to solve both the inviscid (Laplace 
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equation) and viscous flow, whose equations can be written 


aeeording to 


ын gu = = 0 Goan) 
УС - га (5.2) 
where 
р = 1 = constant in laminar flow 
b= 1+ v,/v in turbulent flow 
Four interaction models can be distinguished: Direct, Inverse, 


Semi-inverse, and Simultaneous interaction methods which are 
subject to different boundary conditions. 

The most advanced interaction scheme is the simultaneous 
interaction. We call it the "strong interaction" (direct and 
inverse interactions guarantee weak coupling only). Examples 
in Section V.D are computed bv the Cebeci program using this 
method. Good agreement is obtained between the results of 


interaction schemes and experimental results. 


B. FOUNDATION OF THE INTERACTION SCHEMES 
1. Direct Interaction Scheme 
The direct interaction model is composed of a direct 
inviscid and a direct viscous flow solver (see Figure 5.1а). 


The usual sequence is: 
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Airfoil * 
geometry INVISCID é (xn 


Direcc) 





| 33005 
е: 


VISCOUS 
(Direct) 


* 
6 (57 


a) DIRECT 


182011 
geometry 















АТВ ШЕ 
Direc} 


VISCOUS 
‘Inverse, 





СО? 
Ua OX) Ux) $0? 





C) SEMI-INVERSE 


Ус 6" (x) 
(Inverse) 


ao) 


; 


мо INVISCID 


Ппуегве) 
* 
85 (х) 


COMVERGENCE? 


b) INVERSE 


да трети 
geometry 


Утв 
(Direct) 





' 
| ISTERAGTION ая 





d) SIMULTANEOUS 


Figure 5.1. Global Organization of Interaction Methods 
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(1) Calculate the external velocity distribution by 
inviscid flow computations. 


(2) Calculate the displacement thickness, 6*, by viscous 
flow computations using the external velocity as 
Бойнсагу» Condition.: 


(3) Compute an updated shape of the displacement body and 
repeat steps l and 2 until the results converge. 


Unfortunately, the direct boundary layer method suffers 
numerical breakdown at the point of separation (Goldstein 
singularity). Therefore this scheme is not appropriate to 
handle airfoil flows with separation. 
2. Inverse Interaction Scheme 
This method was introduced to overcome the singularity 
problems near separation. It combines an inverse inviscid and 
an inverse viscous flow solver (see Figure 5.1b). However, 
the overall procedure suffers from very slow convergence. 
Due to this shortcoming one shall apply this inverse scheme 
to regions of separated flow only. 
3. Semi-Inverse Interaction Scheme 
This method combines a direct inviscid flow solver 
with an inverse viscous flow solver with the same input (dis- 
placement pie nese: This leads to two external velocity 


distributions, О ет (СК) апа 0. (x) (see Figure 5.1с). Satis- 


V 


factory convergence is ensured by a relaxation formula, which 


is introduced to define an updated displacement thickness 


СЕКЕ ояз оп . 
Uey UO 
$*ew 5) = бота (Х) (1 T So = 15111 (22157) 
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where ш is a relaxation parameter. The numerical weakness of 

the purely inverse scheme is improved by this method, but 

both inviscid and viscous regions are still coupled loosely. 
4. Simultaneous Interaction Scheme 

The simultaneous interaction scheme emphasizes strong 
interaction between the outer inviscid and the inner viscous 
region. The external velocity U (3) and the displacement 
thickness ó*(x) are treated as unknown quantities. An addi- 
tional relation is added, the so-called interaction law which 
can be given by the "blowing velocity" concept. 

The equations are solved by the inverse method with 
successive sweeps over the airfoil surface (see Figure 5.ld). 
This method is compatible with the weak interaction scheme 
where both inviscid and viscous regions are coupled loosely. 

For each sweep, the external velocity for the boundary 


layer equation is written as 
O 
U(x) = UG X) + SU, (x) (5245) 


where 


02 (х) is the inviscid velocity; 


SU, (x) is the perturbation due to the dis- 
placement effect of a boundary layer. 


The blowing velocity concept is introduced to get the pertur- 


bation velocity SU. by the interaction law. The displacement 
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effect of a boundary layer can be modelled by ejecting fluid 


о ЕО а1:2011'5 surface (see Figure 5.2). 





puo 5.2. Бісоміпа ейссісу СолсерЕ 


A properly arranged source distribution on the surface dis- 
places the streamlines away from the surface such that the 
virtual displacement body becomes a streamline. 

Our first goal is to determine the source strengths such 
that the tangential flow condition on the displacement body 


takes the form 


мае”) 2 
ЖЕЗ dx ERS) 





To achieve this goal we use the thin airfoil approximation: 


(1) The displacement thickness is assumed to be so 
small that u-velocity components do not vary across 
the layer. 


(2) The airfoil in this connection can be represented by 
a straight line. This approximation implies that 
the blowing velocity v (x,0) equals half of the 
Source strength. 
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Therefore, 


ois) = v(x,0) 


ду 
28 За ам пеш 
28/0 ) , ду dy 
ао 
2 аб | е gx 
еба ах 
снн 
= az Ues ) (5.6) 


where (0 6) is defined as blowing velocity. Our second 

goal is to relate the perturbation velocity, OU. to the blowing 
velocity. This process is quite similar to evaluating tangen- 
tial velocities in the panel method. In fact, this is even 
Simpler because of the straight line surface. 


Хү 


D а 964) 

60. = 5 Г. х-Е 9 (5 200 
а 

where the interaction region is limited to a finite range 

ха [Í X < хр. This integral is referred to as the Hilbert 

integral. Rewriting Eq. (5.4), we finally obtain the inter- 


action law. 


E O 1 
0 (х) = 9х) + = Г. шо — (5.8) 


one numerical implementation of Eq. (5.8) requires a discrete 
approximation of the Hilbert integral. This can be performed 
by using the trapezoidal rule. 

The examples in Section V.D demonstrate that this inter- 
action method can give reliable results for flows up to high 
angle of attack, including flows with bubbles and separation. 
C. CONSIDERATION OF BOUNDARY LAYER TRANSITION AND OF 

TURBULENT FLOW MODELLING 

E Transition 

One of the most important parameters to predict the 
пиле lift of an airfoil is the transition point. Boundary 
layer transition is affected by many parameters, for example, 
the pressure distribution (major parameter), the wall roughness 
and the intensity of the free stream turbulence, etc. Because 
of this fact, the theoretical modeling of transition is very 
complicated and one therefore resorts to experimental 
information. 

In the Cebeci program, the following experimental 
correlation formula is used, which was given by Cebeci and 


Smith (1974) as a relation between R, and Ке, dbstransrtion. 


Ө 


R See ee зээг 


х 
tr Xer ЕГ 


22:12) 


where 


Re. = USX/V 
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Вр = О „6/» 


and @ is the momentum thickness. 
2. Turbulent Flow Model 

Unlike laminar flows, turbulent flows have a compli- 
cated time-dependent behavior. It is too difficult to deal 
with the instantaneous properties. Thus, the mean-flow 
properties are applied in turbulent flow. 

The most common mean-flow models are the "eddy- 
viscosity" formula which are based on thin shear layer 


assumptions. 





= райета: = pv, -= (5. 108 


where Уы is related empirically to the mean flow velocity 


gradient and the length scale. In the Cebeci program, y 


is presented by the algebraic eddy-viscosity formulation of 


Cebeci and Smith. 


{0.4у[1-ехр (-у/А)1}7|5 


u 
ши б Ут 
"л> (5.11) 


со 


Л Л (0-4) dy үүл улан 


More detailed descriptions are listed in Reference 7. 


D. EXAMPLES 
The subsequent examples were computed using a program 


developed by Cebeci and coworkers [Ref. 17], on the NPS IBM 370. 
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1. Demonstration of the Program Capabilities 

The velocity profiles on both upper and lower surfaces, 
as well as in the wake, are presented in Figure 5.3. At 
ӘНЕС ава тео attack (aq = 10°), transition occurs very close 
to the leading edge on the upper surface. The boundary layer 
thickness is quite thin in the accelerated flow region (right 
after the leading edge), but it grows thicker farther down- 
stream in the decelerated flow region (near the trailing edge). 
Eventually, we find a small reversed flow region just before 
the trailing edge in this case. The wake zu Оп shows the 
mixing layer which decays with increasing downstream distance. 

Figure 5.4 demonstrates how lift, drag and the loca- 
tion of transition depend on the angle of attack. The skin 
friction drag is dominant at low angles of .attack, the pressure 
ааа ас high angles of attack (see Figure’ 5.4b). 

Figure 5.5 Shows the distributions of the skin friction 
coefficient, displacement thickness and shape parameter in 
dependence of Reynolds number and angle of attack. In the 
attached flow, the skin friction coefficient decreases along 
the downstream direction until the point of transition 
(laminar region), but increases steeply after transition and 
then decreases again because the skin friction is related to 
the slope of the velocity profile, du/dy, at the surface. At 
high angles of attack, transition on the upper surface occurs 
close to the leading edge and the negative skin friction 


coefficient near the trailing edge indicates separated flow. 
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Figure 5.3. Velocity Profiles on the Surface and шэг 
of FX 63-137 Airfoil (a = 10°, Re = 5 *10°) 
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Also the displacement thickness and the shapeparameter increase 
drastically in regions of separated flow (see Figure 5.5a). 

At low Reynolds numbers, laminar flow can separate at 
measchord and reattach (we call this phenomenon a "bubble"). 
Reattachment can often occur if the pressure gradient de- 
creases rapidly soon after separation, so that a strong reversed 
flow is not established. Thus the shear layer reattaches onto 
внеш расе. Accordingly, the displacement thickness andi the 
shape parameter increase in the bubble region (see Figure 
SMOD) | 

2. Comparison with Experimental Results 

The comparison of pressure distributions is shown in 
ВВ. 6. he Overall lift of inviscid calculations devi- 
ates 20% from the experimental results. The interaction method 
improves the accuracy, but the computation still overestimates 
ene lift. The lIift coefficient obtained by Cebeci's program 
1S approximately 10% larger than the measured one (see Reference 
8). 

Figure 5.7 demonstrates that the accuracy of this method 
drops with lower Reynolds numbers. The reasons for this de- 
creased accuracy at low Reynolds numbers are the higher proba- 
bility to separate and the used turbulence model, which seems 
to be inappropriate for low Reynolds numbers. Therefore, low 
Reynolds numbers and high angles of attack pose severe limi- 
tations (see Figure 5.7a & b). The experimental results are 


p ст References / апа 9. 
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On the other hand, experimental measurements should 
not be expected to be exact. Turbulence level and interference 
effects influence the wind tunnel measurements. Figure 5.7c 
shows good agreement between computed and experimental results 
taken from Reference 10, at low Reynolds numbers. 

The location of transition and separation points have 
an important influence on the lift and drag coefficients. 
Figure 5.8 shows that the prediction of laminar separation, 
reattachment, transition and separation points on the airfoil 
surfaces are in reasonable agreement with the experimental 


data taken from Reference 9. 
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УТ. SUMMARY AND RECOMMENDATIONS 


This thesis treats the problem of incompressible two- 
dimensional steady flow past airfoils or airfoil combinations 
at large angles of attack. A panel method was developed to 
compute the inviscid flow over two cylinders, airfoil-flap 
combinations and airfoils in ground effect. In addition, 
Cebeci's viscous/inviscid interaction method was applied to 
several airfoils and compared iwth available experimental 
data. The agreement is generally quite encouraging. The 
calculations show the sensitivity of the computations to 
Reynolds number and transition. More work is required to 
evaluate the potential and limitations of the viscous/inviscid 


interaction method. 
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